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Abstract 

In this paper we present the Burgers equation in its viscous and non-viscous version. After 
submitting, as a motivation, some applications of this paradigmatic equations, we continue 
with the mathematical analysis of them. This will lead us to confront one of the main problems 
linked to non-linear PDE: The appearance of shocks. Finally, the paper concludes presenting 
several numerical schemes for solving these equations and their corresponding implementation 
in MATLAB. 


Contents 


1 Introduction 2 

2 Burgers Model 2 

2.1 Navier Stokes equations simplification. 2 

2.2 Traffic Flow. 3 

3 Mathematical analysis 4 

3.1 Inviscid. 4 

Characteristic method. 4 

Breaking Time. 5 

The Rankine-Hugoniot jump condition. 6 

Note on weak solutions. 6 

Entropy condition. 7 

Examples. 7 

Vanishing viscosity approach. 9 

3.2 Viscid. 9 

The Cole-Hopf transformation. 10 

Heat equation. 10 

Examples. 11 

4 Numerical methods 12 

4.1 Inviscid. 12 

Up-wind nonconservative. 12 

Up-wind conservative. 13 

Lax-Friedrichs. 13 

Lax-Wendroff. 13 

MacCormack. 13 

Godunov. 14 

4.2 Viscid. 14 

Parabolic Method. 14 

4.3 Code. 15 

4.4 Numerical experiments. 18 



























1 INTRODUCTION 


2 


1 Introduction 

In this paper we will consider the viscid Burgers equation to be the nonlinear parabolic pde 

Uf H - nu x — tV'xx (1) 

where e > 0 is the constant of viscosity. This is the simplest pde combining both nonlinear 
propagation effects and diffusive effects. When the right term is removed from (1) we obtain the 
hiperbolic pde 

'U't T x — 0* (2) 

We will refer to (2) as the inviscid Burgers equation. Note that equation (2) can be rewritten in 
the form 

u 2 

u t + [/(•»)]* = 0 with f(u) = (3) 

where is easily recognizable the structure of a scalar hiperbolic conservation law. Many of the 
ideas presented in this paper (relating to mathematical treatment, numerical methods,...) can be 
formulated in the framework of the theory of scalar hiperbolic conservation laws so for general¬ 
ity we will often refer to (2) in the form (3) and the developments will be valid for a general / in (3). 

There is an important connection between the above equations: Equation (2) is the limit as e —>• 0 
of (1). This is true from the formal point of view but also in the most rigorous sense. This fact 
will be studied in more detail in the paragraph called Vanishing viscosity approach. 


2 Burgers Model 

Burgers equations appear often as a simplification of a more complex and sophisticated model. 
Hence it is usually thought as a toy model , namely, a tool that is used to understand some of the 
inside behavior of the general problem. Here we will present two examples. 


2.1 Navier Stokes equations simplification 

Consider the Navier Stokes equations 

/ V-v = 0, (4.1) 

| (pv) t + V • (pvv) + Vp — pV 2 v = 0. (4.2) ' ' 

It is well known that when p is consider to be the density, p the pression, v the velocity and p the 
viscosity of a fluid, equations (4) describe the dynamics of a divergence free (4.1) incompressible 
(pt = 0) flow where gravitational effects are negligible. 


Simplification in (4.2) of the x componet of the velocity vector, which we will call v x , gives 


dv x 
1 ~9t 


pv' 


, dv x dv x 


pv 


Q y X 

dz 


dp 

dx 


p 


d 2 v x 

dx 2 


d 2 v x d 2 v x \ 
dy 2 + dz 2 ) 


If we consider a ID problem with no pressure gradient, the above equation reduces to 


dv x 

"~dt 


pv 


, dv x 
dx 


p 


d 2 v x 

dx 2 


= 0 . 


(5) 


If we use now the traditional variable u rather than v x and take e to be the kinematic viscosity, 
i.e, e = ^, then the last equation becomes just the viscid Burgers equation as it has been presented 
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in (1). 


When the viscosity p of the fluid is almost zero, one could think, as an idealization, to simply 
remove the second-derivative term in (5). This would lead to 


dv x r dv x 

p -» +fw 7h r = ° 


(«) 


which, after making u = v x and dividing by p, becomes the inviscid Burgers equation as it is shown 
in (2). It turns out that, in order to use (6) as a model for the dynamics of an inviscous fluid, it 
has to be supplemented with other physical conditions (section 3.1) which will prevent equation 
(6) from developing physical meaningless solutions. This extension is worth because working with 
(6) is much easier than dealing with (5). 


2.2 Traffic Flow 

Consider the flow of cars on a highway and let p(x,t) denote the density of cars and f(x,t) the 
traffic flow. We will also consider p* to be the restriction of p to a certain range, 0 < p* < pmax, 
where p max is the value at which cars are bumper to bumper. 


Since cars are conserved, the density of cars and the flow must be related by the continuity 
equation 


dp* df ^ 

—— -\ -— = 0 . 

dt* dx* 


(7) 


Obviously, the first expression in which one thinks for the flow is / = vp* where v is the velocity. 
However, it turns out that in order to reflect the fact that drivers will reduce their speed to account 
for an increasing density ahead we should suppose that / is a function of the density gradient as 
well. A simple assumption is to take 


/(p*) = p*v(p*) — where D is a constant. (8) 

We are assuming also that the velocity v is a given function of p*: On a highway we would 
optimally like to drive at some speed v ma x (the speed limit perhaps) but with heavy traffic we slow 
down, with velocity decreasing as density increases. The simplest relation that is aware of this is 

-/>*). (9) 

Pmax 

Putting (8) and (9) into (7) leads to 


dt* dx * 


Vmax 

Pmax 


C p 


max 



= D 


dp* 2 

dx* 2 


( 10 ) 


Scaling through v max = |p = pmaxp *, x = xqx* and t = tot * results in 


pt + [(i - p)p\ x 


£Pxx 


with 


D 

€ = - 

^maxX' 0 


and 


0 < p < 1. 


(ii) 


The transformation u = 2p — 1 leads to the viscid Burgers equation as it is shown in (1) with the 
conditions —1 < u < 1. 
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3 Mathematical analysis 

From the mathematical point of view Burgers equations are a very interesting and suggestive topic: 
It turns out that a study of them leads to most of the ideas that arises in the field of nonlinear 
hiperbolyc waves. 

3.1 Inviscid 

We will focus first on equation (2). Specifically, we will deal with the initial value problem 

f 'U't H - 0? x G K., t 0, /h 

\ u(x, 0) = uo(x), xgM. 

As it as has been suggested previously, although (12) seems to be a very innocent problem a priori 
it hides many unexpected phenomena. 

Characteristic method. The similarity with the advection equation suggests considering, as 
a first approach to solve (12), the characteristic method. In this case the characteristic equation 
would be 


/ x'(t) = u(x(t),t), t > 0, . , 

\ x(0) = x 0 . [ } 

If x(t) and u(x,t) (g C 1 ) are solutions of (13) and (12) respectively, then 

[u(x(t),t)\ = u t (x(t),t) + x'(t)u x (x(t),t) = u t (x(t),t) + u(x(t),t)u x (x(t),i) = 0, 

i.e, u is constant along the characteristic curve x(t) and therefore 

u(x(t),t) = u(x(0), 0) = u 0 (x 0 ), (14) 

which considering the sistem (13), leads to conclude that the characteristic curves are straight lines 
determined by initial data : 

x = xq + uo(xo)t, t > 0. (15) 

In principle, one could invert (15) to obtain xo = xo(x, t). Then, using (14), one would obtain the 
solution u(x,t) = uo(xo(x,t)). However, as the required inversion usually can not be accomplished 
analytically, one could use a symbolic calculation program to construct a discretized solution of 
(12) by dragging the initial data 'Uo(^o) along the characteristic line (15). This strategy is followed 
by the following program in MATHEMATIC A: 


Clear[fO, fval, x, f] 

fO[x_] = Exp[-(2 (x - 1))~2]; (*Condicion inicial*) 
x[t_, x0_] = xO + fO[xO] t; 
f [t_, x0_] = fO[xO] ; 

fval[t_] := Table[{x[t, xO], f[t, xO]}, {xO, -.5, 3, .1}] 

(*Dibujo de las caracteristicas*) 

Plot[Table[x[t, xO], {xO, -.5, 3, .1}], {t, 0, 2},AxesLabel -> {Text[Style ["t", Italic, 23]], 
Text [Style["x", Italic, 23]]}] 

(*Dibujo de la solucion u(x,t)*) 

ListPlot3D[{Table[Table[{x[t, xO], t, f[t, x0]>, {xO, -.5, 3, .1}], {t, 0, 2, 0.1}]}, 
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ColorFunction -> "DarkRainbow", PlotStyle -> Directive[Opacity[0.9]], MeshFunctions -> {#2 &}, 
Mesh -> 5, PlotRange -> All, Axes -> {True, True, True} , 

Boxed -> False, ImageSize -> 600,AxesLabel -> {Text[Style["x", Italic, 23]], 

Text [Style["t", Italic, 23]], Text[Style["u", Italic, 23]]}] 




Figure 1: Output of the programme for uo(x) = e ^ x Some characteristic lines are presented 

in the first picture. 


Looking to this example one quickly finds that problem (12) exibits (under certain initial con¬ 
ditions) what is called the wavebreaking phenomenon: The peak of the pulse moves the fastest 
because wave speed increases with increasing amplitude. This can also be understood by looking 
to the corresponding slope of the characteristic line in the (x, t) plane. Eventually the peak over¬ 
takes the rest of the pulse, or the characteristic cross with another one, and the solution becomes 
multiple valued. The time value ts at which this happens for the first time is called the breaking 
time. 


Breaking Time. In order to determine the breaking time, let us consider two characteristics 
that arise from initial conditions x\ and X 2 = x\ + Ax. According to (15), these characteristics 
will cross when 

X\ + Uo(xi)t = x 2 + Uo(x 2 )t. 


Solving for t leads to 


X\ - X 2 

Uo(x i) - Mo (x 2 ) 


Ax 

Uo(x i) - m 0 (:Xi + Ax)' 


(16) 


When Ax —} 0 then the time in (16) converges to t = — ,} , . To find the breaking time £# we 

u O\ x l) 

find the minimum (positive) value for £, 


ts = m in 



Although the solution obtained via the characteristic method seems to be valid for t < it 
is clear that mathematically multiple valued solutions in a region is unphysical in most of the 
situations (think of u as density) so it can not be accepted for t > £#. What is then the solution 
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of (12) after the breaking time? This question is in the basis of what is called the shock-fitting 
problem: 

The only way to establish a solution after the breaking time is to allow discontinuities of u. Such 
discontinuity is called a shock. This requires some mathematical extension of what we mean by a 
solution to (12), since strictly speaking the derivatives of u will not exist at a discontinuity. It can 
be done through the concept of a weak solution. However, by expanding our class of solutions to 
include discontinuous solutions, we can no longer guarantee the uniqueness of the solution to (12). 
The nonuniqueness can be resolved only by appeal to physically inspired criteria. More specifically, 
we need some shock conditions relating the jumps of u across the discontinuity. 


The Rankine-Hugoniot jump condition. The Rankine-Hugoniot jump condition determines 
the position of a shock at a given time. This condition emerges when one consider the equation 
(3) in integral form by integrating respect to x over X 2 < x < aq, 

IJ u(x,t)dx +f{u) x x \ = 0. (17) 

The physical background is clear, the equation (17) is the traditional form of a conservation law: 
The rate of change of the total amount of u in any section x^ < x < x\ must be balanced by the 
net inflow across x\ and X2 . Suppose now that there is a discontinuity at x = s(t) and that x\ 
and X2 are chosen so that X2 < s(t) < x\. Suppose u and its first derivative are continuous in 
x\ > x > s(t) and in s(t) > x > X2, and have finite limits as x s(t ) from above and below. 
Then, if we assume (17) we have: 


f(u)x 2 - f(u)x 1 




dx = 


U2S — U\S + 



where U 2 ,U\ are the value of u as x s(t ) from below and above respectively and s = 

Since u t is bounded in each of the intervals separately, the integral tend to zero in the limit as 
x\ 5 + ,^2 s~ and, therefore, 


f(u) 2 - /(«) 1 = (u 2 - Ui)s. 

Note that the argument above is valid for a general hiperbolic conservation law of the form 
u t + [f(u)] x = 0. If we focus on equation (3) we obtain 


f(u)2 - f(u) 1 
«2 - «1 



M 2 - Ml 


-(Ml +M 2 ) 


(18) 


This is the first condition that the shock of our discontinuous single valued solution will have to 
satisfy. 


Note on weak solutions. In this section we will work within hiperbolic conservation laws, i.e., 
as it has been told in the introduction, equations in the form 


Ut + [f(u)] x = 0 


( 19 ) 
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with some / (see (3)). Mathematically, a discontinuous solution with continuously differentiable 
parts satisfying (19) and such that s = can be considered a weak solution of (19). As 

a motivation for the idea, let </>(x, t) be an infinitely smooth function in M x [0, +oc) that vanishes 
on the boundary. Let u(x, t ) be a smooth solution of the hyperbolic conservation law given. Then 

o= r 

J 0 l-oo 


> = 0 on bound. 


/ OO p OO pOO pOO 

(<H~ dx+ WWl-oodt- / / {4> t u + dxdt 

-oo JO Jo J — oo 

pOO poo 

- / (4>tU + 4>xf(u)) dxdt 

J 0 J — oo 

In view of this development it is natural to say that u(x,t) is a weak solution of the conservation 
law if for any 0(x, t) the equation 



{(j) t u + 4> x f(u )) dxdt = 0 


( 20 ) 


holds. If we consider now a weak solution of (19) with a simple jump discontinuity across 
Q = {(s(t),t)} it is easy to show, dividing the integral domain in (20) by f2, that s = . 


At first sight the weak solution concept appears to be the appropriate mathematical tool to work 
with. However, it turns out that there are still some situations in which this concept is not enough 
to guarantee the uniqueness. Becouse of that, there is an additional condition, the so called entropy 
condition, introduced to eliminate nonphysical weak solutions. There is a number of variations in 
which this condition can be presented, we will mention only the simplest one. 

Entropy condition. For the equation (19), a discontinuity propagating with speed s = 
satisfies the entropy condition if f{u 2 ) > s > f'{ui). For the Burgers equation this entropy con¬ 
dition reduces to the requirement that if a discontinuity is propagating with speed s then u 2 > U\. 

As its name suggests, this condition is the mathematical translation of the condition that says 
that in every physical process the entropy of the system is nondecreasing. As is well known, this 
is a fundamental assumption in thermodynamics. 


Examples. 


1. Riemann problem. The initial value problem (12) with discontinuous initial condition of 
the form 


U 0 (x) 


u l , X < 0 , 

Ur, X>0 


is called a Riemann problem. As a first example we will take the simple assumption ul = 1 
and ur = 0. If we use the characteristic method presented above we obtain 
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Figure 2: 


It is obvious that this solution needs a shock fitting just from the beginnig. Obeying the 
Rankine-Hugoniot condition we obtain that the discontinuity must travel with the speed 
s = \ and therefore the solution becomes 


u 



1 , 

0 , 


X — 2 ’ 
X > 2 * 


Note that this solution also satisfies the entropy condition. In this case it can be shown that, 
in fact, this is the unique weak solution for the problem. 


If we consider now just the opposite problem, i.e, Ur = 0 and Ur = 1, the hole situation 
changes. There ere at least two solutions satisfying the Rankine-Hugoniot condition: 




0, x < 0, 

|, 0 < x < £, and U 2 (x,t) 

1, x > t. 


0, X < f, 
1, x > §. 


In this situation, the fact that the second one does not satisfy the entropy condition allows 
us to discard it and keep the continuous solution u\. 


2. In order to get a more accurate idea of what it means to fit a shock we will consider the 
initial condition 

{ 1, x < 0 

1 — x, 0 < x < 1 (21) 

0, x > 1 

for (12). If we draw some characteristics lines we obtain 



Figure 3: 
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As it is seen in the figure we have £# = 1 so it is from there that we have to fix the solution. 
Appealing again to the R-H condition, the shock path for t > ts must be s(t) = In 

the following figure it is shown the solution to the problem: Until ts = 1 we can keep the 
solution obtained via the characteristic method but then we have to rid it and fix it. 



Figure 4: In this figure we can see the part of the characteristic solution that has to be replaced 
in order to obtain the solution we are lookin for (the darker part of the figure). 


Note that the entropy condition is satisfied by this solution. 

Vanishing viscosity approach. Fitting a shock in the region of overlapped characteristics fol¬ 
lowing the Rankine-Hugoniot and the entropy condition seems mathematical unnatural. However, 
there is a more plausible and natural approach to construct the discontinuous entropy solution. 
This is the so called vanishing viscosity approach. The first step is to add the viscosity-dispersion 
term eu xx to obtain equation (1). In the next section we will see that adding this term has an 
important effect: it suppresses the wavebreaking. The basic reason is easy to understand: disper¬ 
sion causes waves to spread and this acts against the steepening effect of the nonlinearity. Hence 
we expect to obtain a smooth solution of (1) even with discontinuous initial data. The vanishing 
viscosity approach is based on the natural intuition that this smooth solution of (1) will approach 
a shock wave as e —>> 0. On the physical interpretation point of view this is the natural way of 
setting the correct solution of problem (12): As it has been said in the subsection (2.1), the inviscid 
Burgers equation usually appears as an idealization of a more precise model where there is always 
some degree of viscosity. That is why the only relevant solutions are those that can be obtained 
with this approach. There is a result due to Kruzkov that guarantees the uniqueness of the solution 
obtained by this way. 

The vanishing viscosity approach would be completely useless if one can not find the solution to 
the correspondig initial value problem (the problem (12) where it has been added the viscosity- 
dispersion term). Fortunately the Cole-Hopf transformation, introduced by Cole (1951) and Hopf 
(1950) independently, transforms the equation (1) into the heat equation and therefore it can be 
solved. The Cole-Hopf transformation is introduced in the next section. 

3.2 Viscid 

Consider now the initial value problem for the viscid Burgers equation 

j u t + uu x = eu xx , xGM, t > 0, e > 0, 

| u(x,0) = uq(x), xgM. 


(22) 
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The Cole-Hopf transformation. The Cole-Hopf transformation is defined by 

0 Tx 
u = —2e —. 


Operating in (23) we find that 


2e (tpttpx - <y 0 (p xt 


4 e 2 ip x (w ra - ifil 


dUxX - 


2e [2p x 3(p(f> X xTx T T Txxx) 


Substituting this expressions into (1), 

2e ( ppxt T Tx {}Pt ^Txx) T tppxx 


— 0 4 > ppxt T Tx i}Pt ^Txx) T tppxxx — 0 


(fx (ipt - ZVxx) = <p (<Pxt - eipxxx) = V(Vt~ Z<Pxx)x • 

Therefore, if p solves the heat equation p t — dp xx = 0 x E M, then u(x, t ) given by transformation 
(23) solves the viscid Burgers equation (1). 

To completely transform the problem (22) we still have to work with the initial condition function. 
To do this, note that (23) can be written as 


u = —2d (log p) x ' 


Hence 


It is clear from (23) that multiplying p by a constant does not affect u, so we can write the last 
equation as 

= ^ } dy ). (24) 

The initial condition on (22) must therefore be transformed by using (24) into 

(_ rx u 0 (y) , \ 

p(x, 0) = po(x) — 0 2e V ' • 

In summary, we have reduced the problem (22) to this one 

f p t - £Txx =0, x G M, t > 0, e > 0, 

\ (_ rx upjy) i \ ( 25 ) 

1 p(x, 0) = po(x) = ev 2e y ) . x G M. 


Heat equation. The general solution of the initial value problem for the heat equation is well 
known and can be handled by a variety of methods. Taking the Fourier transform with respect to 
x for both heat equation and the initial condition po(x) in problem (25) we obtain the first order 


( Tt = £ 2 e<^, £ E M t > 0 e>0 

l £(£,0) = £ 0 (O l e R, 

where ip(£,t ) = f^° ip(x,t)e l ^ x dx. The solution for this problem is 


<p(£,t) = (po{d)e^ et . 
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To recover (p(x,t) we have to use the inverse Fourier transformation F 1 , namely, 

</>0 ,t) = = F _1 (^ 0 (O ei2et ) = <Po(x) * F~ 1 {e^ tt ), 

where * denotes the convolution product. 

On the other hand 

_ -i , £-2 j 1 x ^ 

F 1 (e* et ) = —^=e 4 et 
V ' 2 V^t 

so the initial value problem (25) has the analytic solution 

1 

M, t ) = —r= MO 

2\/7ret J -oo 


e 4et at 


Finally, from (23), we obtain the analytic solution for the problem ( 22 ) 


u(x, t ) 


rl l 


(27) 


Examples. We will use now (27) to draw the exact solution of (22) with different initial condi¬ 
tions. In particular we will take the same initial conditions that we took for the inviscid case and 
we will vary the viscosity parameter e in order to visualize the vanishing viscosity effect. 




Figure 5: uq(x) = e ^ x with e = 0.1 and e = 0.01 respectively. 


2 . 



Figure 6 : Riemann problem with ul = 1, ur = 0 and e = 0.1 and e = 0.01 respectively. 
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Figure 7: Same initial condition that the second example of last section with e = 0.1 and e = 0.01 
respectively. We can see how the solution approaches the fixed one (figure 4) as e -A 0. 


4 Numerical methods 

At this point we should be convinced of the complexity that nonlinearity hides from the point 
of view of the mathematical analysis. This complexity also arises when one trie to solve Burgers 
equation using numerical methods. Major problems arise when one tries to approximate the 
solutions for which we have admit discontinuities: It could happen that the method converges to 
another weak solution of our original equation (or that is the wrong weak solution, i.e, does not 
satisfy the entropy condition). In order to avoid that to happend, numerical schemes must satisfy 
certain non-obvious conditions. 

4.1 Inviscid 

Up-wind nonconservative. If we consider the inviscid burgers equation in the quasilinear 
form 

u t + uu x = 0, (28) 

then a natural finite difference method obtained by a forward in time and backward in space 
discretization of the derivatives is 


jjn+l = jjn _ 




u?- 1 )' 


(29) 


We will refer to (29) as the Up-wind nonconservative scheme. Although this method is consistent 
with (4.1) and is adequate for smooth solutions, it will not converge in general to a discontinuous 
weak solution as the grid is refine. 


To prevent a method from converging to non-solutions, there is a simple condition that we can 
require : The method will have to be in conservation form, i.e, in the form 


ur 1 


V? 


U l F U?-P+1> ->U? +q )~ F ( U ?-P- 1 . u ?~ 


P’ • 


VjVi)] 


(30) 


where F is some function of p + q + 1 arguments called the numerical flux function. Methods 
that conform to this scheme are called conservative methods. The following methods belong all to 
this class. We will first present the schemes for a general scalar conservation law u t + [f(u)\ x = 0 
and then particularize them to the inviscid Burgers equation showing how thay are going to be 
implemented in matlab. 
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Up-wind conservative. If we consider a general scalar conservation law 

u t + [/(«)]* = 0 

and use the standard finite difference discretizations, we obtain the conservation method called 
Up-wind conservative, namely, 

Vl H = V ?-t[f<Vj)-f < 31 ) 

Note that this is in the form (30) with p = 0, q = 1 and F(U, V) = f ( U ). 


For the Burgers equation we have 


C? +1 = v? -1 \l Kf 



(32) 


Lax-Friedrichs. The Lax-Friedrichs method to nonlinear systems takes the form 

f(r w = \vu + k +1 ) - a [/(t? +I ) - f&?_,)] 

Note that this is in the form (30) with p = 0, q = 1 and F(U, V) = ^(U — V) + \ {f{U) + f(V)). 


For the Burgers equation we have 


n-f- l _ 






5 (^.) J - 5 (^.) 2 


(33) 


Lax-Wendroff. The methods consider hitherto are all of first order. The Lax-Wendroff method 
to nonlinear conservation laws is a second order method and takes the form 




TT n _ 
u j+l 


^(/W + .)-/W- i))+S 


2 h 2 L 


Vi (/(%) - /(^)) - Vi CW) - 


where A.+ i is the jacobian matrix A(u) = f(u) evaluated at + LT+ 1 ). 

For Burgers equation we have /'(iz) = u so 



(34) 


MacCormack. Another method of the same type is known as MacCormack’s method. This 
method uses first forward differencing and then backward differencing to achieve second order 
accuracy: 


u* = U?A[f(U? + 1 )-f(Up] 

u ? +1 = \v? + up-L[f(uD-f{up i)] 


( 35 ) 
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Godunov. The last method for solving Burgers equation that will be presented in this paper 
belongs to the so called finite volume methods. The idea of Godunov’s method is the following. Let 
Uj 1 be a numerical solution on the n-th layer. Then we define a function u n (x , t) for t n < t < t n +1 
as follows. At t = t n , 


u n (x,t) mu?, Xj - ^ < x < Xj + j = 2,..., n - 1. 

Then we define u(x, t ) to be the solution of the collection of Riemann problems on the interval 
[t n ,t n + 1 \- If k is small enough so that the characteristics starting at the points Xjt_ | do not inter¬ 
sect within this interval, then u(x, t n + 1 ) is determined unambiguously. Then the numerical solution 
on the next layer, U ™ +1 is defined by averaging u(x, t n+ i) over the intervals Xj — | < x < Xj + |. 
This idea reduces to a very simple conservative scheme: 



+ 

II 

1 

£-1 

[F(U?,U? +1 )-F u?)} 

/ * \2 

where the numerical flux F is defined by F(U,V) = where iz* is defined as follows 

If U > V then 

«* = j 

f U, if u + v > 0, 

[ V, in other case. 

If U < V then 

* J 

u = < 

( U, if C/ > 0, 

V, if V < 0, 

0, if t/ < 0 < V. 



V 


4.2 Viscid 

Parabolic Method. Consider equation (1) in the form 

u 2 

u tT[f(u)] x = eu xx with f(u) = 

If we formally integrate (36) from to Xj + i and rewrite the equation, we obtain 


X J + 1/2 


X 3 — 1 / 2 


ut dx - [eu x ]*'+J£ 


[/(«)] 


x j + 1/2 
X 3 — 1/2 


(36) 


(37) 


We will now look for approximations of each of the terms in (37). We have that 

r x j+ 1/2 d u 

/ u t dx &—(xj,t)h, 

4 x 3- 1/2 


dt 


- i eu ^Z + -i/l = e K( x i- 1 / 2 .*) - M*(^'+i/ 2 ,i)] ~ e 


7i(xj , t) — iz(#j_i, t) , t) — u(Xj , t) 

h h 


u(xj+ 1 , t) — 2 u{xj,t) + u(xj- i,t) 
h 


and 


- :/(»)];^i;(; = /(»(•<•.,• 1/2-0) - /<"(•'•.,>1/2-0) 
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Substituting these expressions into (37), dividing by h and taking Uj(t ) to be the function u(xj,t), 
we obtain the system of ordinary differential equations 

dUj U j+ 1 - 2 Uj + Uj -1 = f( u j- 1) - f(U j+ 1) 

dt £ h 2 ~ h 

As a computational method we discretize the time derivative in the last expression by a forward 
difference to obtain the explicit method 


Un+i = u? + k e 


U ? + 1 - 2 Uf +UJ_WJ+i) ~ fU-'j ,)' 


h 2 




(38) 


where we take f(U-+ 1 ) to be the average of f{Uj) and /( U, + -|)■ 


4.3 Code 


°/ 0 -Burgers program- 

°/ 0 This program computes numerical solutions for viscid and inviscid Burgers equation. 

°/ 0 It needs the functions df.m, f.m, nf.m and uinit.m . 

clear all; 

°/ 0 Selection of equation and method. 
type_burger = menu( 5 Choose the equation: 5 , ... 

5 Inviscid Burgers equation 5 ,’Viscid Burgers equation 5 ); 

if (type_burger == 1) 

method = menu(’Choose a numerical method: 5 , ... 

5 Up-wind nonconservative 5 , 5 Up-wind conservative 5 ,’Lax-Friedrichs 5 , 5 Lax-Wendroff 5 ,’MacCormack 5 ,’Godunov 5 ); 
else 

method = menu(’Choose a numerical method: 5 , ... 

’Parabolic Method 5 ); 
end 

ictype = menu(’Choose the initial condition type: 5 , ... 

’Piecewise constant (shock) 5 ,’Piecewise constant (expansion) 5 , 5 Gaussian 5 , 5 Piecewise continuous 5 ); 

“/Selection of numerical parameters (time step, grid spacing, etc...), 
xend = 2; °/ 0 x-axis size, 

tend = 2; °/ 0 t-axis size. 

N = input(’Enter the number of grid points: ’); 

dx = xend/N; °/ 0 Grid spacing 

dt = input(’Enter time step dt: ’); 

x = 0:dx:xend; 
nt = floor(tend/dt); 
dt = tend / nt; 

°/ 0 Set up the initial solution values. 

uO = uinit (x, ictype) ; °/ 0 Call to the function "uinit". 
u = uO; 
unew = 0*u; 

“/Implementation of the numerical methods, 
if (type_burger == 1) 
for i = 1 : nt, 
switch method 
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case 1 °/ 0 Up-wind nonconservative 

unew(2:end) = u(2:end) - dt/dx * u(2:end) . * (u(2:end) - u(l:end-1)); 
unew(l) = u(l); °/ 0 u(3:end): subvector de u desde 3 hasta el final. 

case 2 °/ 0 Up-wind conservstive 

unew(2:end)= u(2:end) - dt/dx * (f(u(2:end)) - f(u(l:end-1))); 
unew(l) = u(l); 

case 3 °/ 0 Lax-Friedrichs 

unew(2:end-1) = 0.5*(u(3:end)+u(l:end-2)) - 0.5*dt/dx * ... 

(f(u(3:end)) - f(u(l:end-2))); 
unew(l) = u(l); 
unew(end) = u(end); 

case 4 °/ 0 Lax-Wendroff 
unew(2:end-1) = u(2:end-l) ... 

- 0.5*dt/dx * (f(u(3:end)) - f(u(l:end-2))) ... 

+ 0.5*(dt/dx)~2 * ... 

( df(0.5*(u(3:end) + u(2:end-l))) .* (f(u(3:end)) - f(u(2:end-1))) - ... 
df(0.5*(u(2:end-1) + u(l:end-2))) .* (f(u(2:end-1)) - f(u(l:end-2))) ); 
unew(l) = u(l); 
unew(end) = u(end); 

case 5 °/ 0 MacCormack 

us = u(l:end-1) - dt/dx * (f(u(2:end)) - f(u(l:end-1))); 
unew(2:end-l)= 0.5*(u(2:end-1) + us(2:end)) - ... 

0.5*dt/dx* (f(us(2:end)) - f(us(1:end-1))); 
unew(l) = u(l); 
unew(end) = u(end); 

case 6 °/ 0 Godunov 

unew(2:end-1) =u(2:end-l)- dt/dx*(nf(u(2:end-1),u(3:end)) - nf(u(l:end-2),u(2:end-1))); 
unew(l) = u(l); 
unew(end) = u(end); 

end 

u = unew; 

U(i,:) = u(:); 
end 

else 

for i = 1 : nt, 
switch method 

case 1 °/ 0 Parabolic method 
D=0.01; 

fminus = 0.5*( f(u(2:end-1)) + f(u(l:end-2)) ); 
fplus = 0.5*( f(u(2:end-1)) + f(u(3:end)) ); 

unew(2:end-1) = u(2:end-l) + dt*(D*(u(3:end)-2*u(2:end-l)+u(l:end-2))/(dx)~2 ... 

- (fplus-fminus)/dx ); 
unew(l) = u(l); 
unew(end) = u(end); 

end 

u = unew; 

U(i,:) = u(:); 
end 


end 
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U=[uO;U]; 

T=0:dt:tend; 

•/.Plot of the solutions, 
figure(1) 
surf(x,T,U) 
shading interp 

xlabel( 5 x 5 ), ylabel( 5 t 5 ), zlabel ( 5 u(x,t) 5 ); 
grid on 

colormap( 5 Gray 5 ); 


The initial condition is introduced through uinit. m: 


function ui = uinit( x, ictype ) 
xshift = 0; 

if (ictype==l) °/ 0 Shock (uL > uR) 

uL = 1; 

uR = 0; 

ui = uR + (uL-uR) * ((x-xshift) <= 0.0); 

elseif (ictype==2) °/ 0 Expansion (uL < uR) 
uL = 0.5; 

uR = 1.0; 

ui = uR + (uL-uR) * ((x-xshift) <= 0.0); 

elseif (ictype==3) °/ 0 Gaussian 
ui = exp(-2*(x - 1)."2); 

elseif (ictype==4) “/Piecewise continuous 
for i=l:length(x) 
if (x(i) < 0) 
ui(i)=1; 

elseif (0<= x(i) && x(i) < 1) 
ui(i)=l-x(i); 
elseif (x(i) >= 1) 
ui(i)=0; 

end 

end 

end 


The functions f .m, df .m and nf .m (used on Godunov): 


function ret = f( u ) 
ret = 0.5 * u."2; 

function ret = df( u ) 
ret = u; 

function ret = nf( u, v ) 
for i = l:length(u) 
if (u(i) >= v(i)) 

if ((u(i)+v(i))/2 > 0) 
ustar(i)=u(i); 
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else 

ustar(i)=v(i); 

end 

else 

if (u(i)>0) 

ustar(i)=u(i); 
elseif (v(i)<0) 

ustar(i)=v(i); 

else 

ustar(i)=0; 

end 

end 

end 

ret =f(ustar); 


4.4 Numerical experiments 

1. Riemann problem with ul = 1 and ur = 0. 





Figure 8: Up-wind non conservative, Up-wind conservative, Lax-Friedrichs, Lax-Wendroff, Mac- 
Cormack and Godunov respectively. 
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2. The inviscid Burgers equation with the initial condition (21) solved by Godunov method: 



Figure 9: 


Note that the obtained solution agrees with the one presented in the second example of 
section 3. 

3. The next examples have been obtained using the parabolic method with N = 100 and 
dt = 0.001. In order to display the effectiveness of the method we can compare the results 
with the exact solutions shown in section 3.2. 



Figure 10: uq(x ) = e ^ x with e = 0.1 and e = 0.01 respectively. 
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Figure 11: Riemann problem with ur — 1, ur = 0 and e = 0.1 and e = 0.01 respectively. 



Figure 12: Same initial condition that the second example of last section with e = 0.1 and e = 0.01 
respectively. 
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